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ABSTRACT 

We present two-dimensional inviscid hydrodynamic simulations of overstable inertial- 
acoustic oscillation modes (p-modes) in black-hole accretion discs. These global spiral 
waves are trapped in the inner- most region of the disc, and are driven overstable by 
wave absorption at the corotation resonance {tc) when the gradient of the background 
disc vortensity (vorticity divided by surface density) at Vc is positive and the disc 
inner boundary is sufficiently reflective. Previous linear calculations have shown that 
the growth rates of these modes can be as high as 10% of the rotation frequency at 
the disc inner edge. We conffim these linear growth rates and the primary disc oscil- 
lation frequencies in our simulations when the mode amplitude undergoes exponential 
growth. We show that the mode growth saturates when the radial velocity perturba- 
tion becomes comparable to the disc sound speed. During the saturation stage, the 
primary disc oscillation frequency differs only slightly (by less than a few percent) from 
the linear mode frequency. Sharp features in the fluid velocity profiles at this stage 
suggest that the saturation results from nonlinear wave steepening and mode-mode 
interactions. 
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1 INTRODUCTION 

In several recent papers (Lai & Tsang 2009; Tsang & Lai 
2009c; Fu k Lai 2011; Horak & Lai 2013), we have presented 
detailed study of the linear instability of non-axisymmetric 
inertial-acoustic modes [also called p-modes; see Kato (2001) 
and Wagoner (2008) for review] trapped in the inner-most 
region of black-hole (BH) accretion discs. This global insta- 
bility arises because of wave absorption at the corotation res- 
onance (where the wave pattern rotation frequency matches 
the background disc rotation rate) and requires that the disc 
vortensity has a positive gradient at the corotation radius 
(see Narayan et al. 1987, Tsang & Lai 2008 and references 
therein). The disc vortensity (vorticity divided by surface 
density) is given by 



c = 



(1) 



where Q.{r) is the disc rotation frequency, ^{r) is the radial 
epicyclic frequency and S(r) is the surface density 0- Gen- 
eral relativist ic (GR) effect plays an important role in the 
instability: For a Newtonian disc, with — (x. r~^^'^ and 
relatively flat 5](r) profile, we have d(^/dr < 0, so the corota- 
tional wave absorption leads to mode damping. By contrast, 
K is non-monotonic near a BH (e.g., for a Schwarzschild BH, 
K reaches a maximum at r = 8GM/c^ and goes to zero at 
nsco = GCM/c^), the vortensity is also non-monotonic. 
Thus, p-modes with frequencies such that d(^/dr > at 
the corotation resonance are overstable. Our calculations 
based on several disc models and Paczynski-Witta pseudo- 
Newtonian potential (Lai & Tsang 2009, Tsang & Lai 2009c) 
and full GR (Horak & Lai 2013) showed that the lowest- 
order p-modes with m = 2, 3, 4, • • • have the largest growth 
rates, with the mode frequencies u ^ /3mQisco (thus giv- 
ing commensurate frequency ratio 2 : 3 : 4, • • • ), where the 
dimensionless constant (3^1 depends weakly the disc prop- 
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^ Equation ([T} applies to barotropic discs in Newtonian (or 
pseudo-Newtonina) theory. See Tsang &; Lai (2009c) for non- 
barotropic discs and Horak & Lai (2013) for full general rela- 
tivistic expression. 
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erties. These overstable p-modes could potentially explain 
the High-frequency Quasi-Periodic Oscillations (HFPQOs) 
observed in BH X-ray binaries (e.g., Remillard & McClin- 
tock 2006; Belloni et al. 2012). 

The effects of magnetic fields on the oscillation modes 
of BH accretion discs have been investigated by Fu & Lai 
(2009, 2011, 2012) and Yu & Lai (2013). Fu & Lai (2009) 
showed that the the basic wave properties (e.g., propagation 
diagram) of p-modes are not strongly affected by disc mag- 
netic fields, and it is likely that these p-modes are robust in 
the presence of disc turbulence (see Arras, Blaes & Turner 
2006; Reynolds & Miller 2009). By contrast, other diskoseis- 
mic modes with vertical structure (such as g-modes and c- 
modes) may be easily "destroyed" by the magnetic field (Fu 
& Lai 2009) or suffer damping due to corotation resonance 
(Kato 2003; Li et al. 2003; Tsang & Lai 2009a). Although 
a modest toroidal disc magnetic field tends to reduce the 
growth rate of the p-mode (Fu & Lai 2011), a large-scale 
poloidal field can enhance the instability (Yu & Lai 2013; 
see Tagger & Pallet 1999; Tagger & Varniere 2006). The p- 
modes are also influenced by the magnetosphere that may 
exist inside the disc inner edge (Fu & Lai 2012). 

So far our published works are based on linear analysis. 
While these are useful for identifying the key physics and 
issues, the nonlinear evolution and saturation of the mode 
growth can only be studied by numerical simulations. It is 
known that fluid perturbations near the corotation resonane 
are particularly prone to become nonlinear (e.g., Balmforth 
& Korycansky 2001; Ogilvie & Lubow 2003). Moreover, real 
accretion discs are more complex than any semi-analytic 
models considered in our previous works. Numerical MHD 
simulations (including GRMHD) are playing an increasingly 
important role in unraveling the nature of BH accretion 
flows (e.g., De Villiers & Hawley 2003; Machida & Mat- 
sumoto 2003; Fragile et al. 2007; Noble et al. 2009, 2011; 
Reynolds & Miller 2009; Beckwith et al. 2008, 2009; Mosci- 
brodzka et al. 2009; Penna et al. 2010; Kulkarni et al. 2011; 
Hawley et al. 2011; O'Neill et al. 2011; Dolence et al. 2012; 
McKinney et al. 2012; Henisey et al. 2012). Despite much 
progress, global GRMHD simulations still lag far behind ob- 
servations, and so far they have not revealed clear signs of 
HFQPOs that are directly comparable with the observations 
of BH X-ray binaries 0. If the corotation instability and its 
magnetic counterparts studied in our recent papers play a 
role in HFQPOs, the length-scale involved would be small 
and a proper treatment of flow boundary conditions is im- 
portant. It is necessary to carry out "controlled" numerical 
experiments to capture and evaluate these subtle effects. 

In this paper, we use two-dimensional hydrodynamic 



^ Henisey et al. (2009, 2012) found evidence of excitation of wave 
modes in simulations of tilted BH accretion disks. Hydrodynamic 
simulations using Q;-viscosity (Chan 2009; O'Neill et al. 2009) 
showed wave generation in the inner disk region by viscous insta- 
bility (Kato 2001). The MHD simulations by O'Neill et al. (2011) 
revealed possible LFQPOs due to disk dynamo cycles. Dolence 
et al. (2012) reported transient QPOs in the numerical models 
of radiatively inefiicent flows for Sgr A"*". McKinney et al. (2012) 
found QPO signatures associated with the interface between the 
disc inflow and the bulging jet magnetosphere (see Fu &; Lai 2012). 
Note that the observed HFQPOs are much weaker than LFQPOs, 
therefore much more difficult to obtain by brute- force simulations. 



simulations to investigate the nonlinear evolution of coro- 
tational instability of p-modes. Our 2D model has obvious 
limitations. For example it does not include disc magnetic 
field and turbulence. However, we emphasize that since the 
p-modes we are studying are 2D density waves with no ver- 
tical structure, their basic radial "shapes" and real frequen- 
cies may be qualitatively unaffected by the turbulence (see 
Arras et al. 2006; Reynolds & Miller 2009; Fu & Lai 2009). 
Indeed, several local simulations have indicated that den- 
sity waves can propagate in the presence of MRI turbulence 
(Gardiner & Stone 2005; Fromang et al. 2007; Heinemann 
& Papaloizou 2009). Our goal here is to investigate the sat- 
uration of overstable p-modes and the their nonlinear be- 
haviours. 



2 NUMERICAL SETUP 

Our accretion disc is assumed to be inviscid and geometri- 
cally thin so that the hydrdynamical equations can be re- 
duced to two-dimension with vertically integrated quanti- 
ties. We adopt an isothermal equation of state throughout 
this study, i.e. P — c^S where P is the vertically integrated 
pressure, S is the surface density and Cs is the constant 
sound speed. Self- gravity and magnetic fields are neglected. 

We use the Paczynski-Witta Pseudo-Newtonian poten- 
tial (Paczynski & Witta 1980) to mimic the GR effect: 



GM 
r — rs^ 



(2) 



where rs = 2GM/c^ is the Schwarzschild radius. The corre- 
sponding Keplerian rotation frequency and radial epicyclic 
frequency are 



GM 1 



■ rs 



- 3rs 



- rs 



(3) 
(4) 



In our compuation, we will adopt the units such that the 
inner disc radius (at the Inner-most Stable Circular Orbit 
or ISCO) is at r = 1.0 and the Keperian frequency at the 
ISCO is ^^isco = 1. In these units, Ts = 1/3, and 



3r 



(5) 



Our computation domain extends from r = 1.0 to r = 4.0 
in the radial direction and from = to = 27r in the 
azimuthal direction. We also use the Keplerian orbital period 
(T = 27t/Qk = 27r) at r = 1 as the unit for time. The 
equilibrium state of the disk is axisymmetric. The surface 
density profile has a simple power-law form 



(6) 



which leads to a positive vortensity gradient in the inner 
disc region. The equilibrium rotation frequency of the disc 
is given by 



^^o(r) 



4/9 



r(r - 1/3)2 



(7) 



Throughout our simulation, we will adopt Cg — 0.1 so that 
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Figure 1. Evolution of the radial velocity amplitude |wr|max (evaluated at r = 1.1) for three runs with initial azimuthal mode number 
m = 2 (left panel), m = 3 (middle panel) and with random perturbations (right panel). The dashed lines are the fits for the exponential 
growth stage (between ~ 10 and ~ 30 orbits) of the mode amplitude. 




Figure 2. Comparison of the radial profiles of velocity perturbations from non-linear simulation and linear mode calculation. The top 
and bottom panels show the radial and azimuthal velocity perturbations, respectively. The left and right panels are for cases with 
azimuthal mode number m = 2 and m = 3, respectively. In each panel, the dashed line is taken from the real part of the complex 
wavefunction obtained in linear mode calculation, while the solid line is from the non-linear simulation during the exponential growth 
stage (at T = 20 orbits), with the quantities evaluated at (j) = OAtt. Note that the normalization factor is given in the y-axis label for 
the nonlinear simulation results. 



We solve the Euler equations that govern the dynam- 
ics of the disc flow with the PLUTO code (Mignone et 

^ publicly available at http: / /plutocode.ph.unito.it/l 
© 2012 RAS, MNRAS 000, [l]-?? 



al. 2007), which is a Godunov-type code with multiphysics 
and multialgorithm modules. For this study, we choose a 
Runge-Kutta scheme (for time integration) and piecewise 
linear reconstruction (for space integration) to achieve sec- 
ond order accuracy, and Roe solver as the solution of Rie- 
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mann problems. The gird resolution we adopt is {Nr x A^^) = 
(1024 X 2048) so that each grid cell is almost a square. Our 
runs typically last for 100 orbits (Keplerian orbits at inner 
disc boundary). To compare with the linear mode calcula- 
tions (Lai & Tsang 2009), the inner disk boundary is set 
to be reflective with zero radial velocity. At the outer disc 
boundary, we adopt the non- reflective boundary condition. 
This is realized by employing the wave damping method (de 
Val-Borro et al. 2006) to reduce wave reflection. This im- 
plementation mimics the outgoing radiative boundary con- 
dition used in the linear analysis. 



3 RESULTS 

We carried out simulations with two different types of initial 
conditions for the surface density perturbation. In the first, 
we choose ^S(r, cj)) to be randomly distributed in r and (j). In 
the second, we impose ^S(r, 0) oc cos(m(/)) that is randomly 
distributed in r, so that the perturbation has an azimuthal 
number m. In all cases, the initial surface density perturba- 
tion has a small amplitude (I^S/Sol ^ 10~^). 

Fig. [1] shows the evolution of the radial velocity ampli- 
tude near the inner disc radius for runs with initial azimuthal 
mode number m = 2, m = 3 and random initial perturba- 
tion, respectively. This velocity amplitude is obtained by 
searching for the maximum \ur\ at a given r by varying (j) 
for each given time point. We chose r = 1.1 because this 
is where the largest radial velocity perturbation is located 
(see the upper panels of Fig. [2]). We see that in all cases 
there are three stages in the amplitude evolution. The first 
stage occupies roughly the first 10 orbits, during which the 
initial perturbation starts to affect the flow and presumably 
excites many modes/oscillations in the disc. In the second 
stage (from ^ 10 orbits to ~ 30 orbits), the fastest growing 
mode becomes dominant and undergoes exponential growth 
with its amplitude increased by about four orders of mag- 
nitude. In the last stage (beyond ^ 30 orbits), the pertur- 
bation growth saturates and its amplitude remains at ap- 
proximately the same level. A fit to the exponential growth 
stage gives the growth rate of the fastest growing perturba- 
tion, which is the slope of the fitted straight line. For the 
m — 2 run, we find 0.17/ log ]^Q(e)/27r 0.0637 (in units of 
the orbital frequency at r = 1) as the growth rate, which is 
quite consistent with the result from our linear eigenmode 
calculation uji ^ 0.0632 (the imaginary part of the eigen- 
frequency; Lai & Tsang 2009; Fu & Lai 2011; see Table [J). 
For the m = 3 run, our simulation gives 0.074 as the mode 
growth rate, close to Ui = 0.0733 from the linear calculation. 
In Fig. [21 we plot the radial profile of the velocity perturba- 
tion at a given time during the exponential growth stage of 
the simulation, and we compare it with the wavefunctions 
obtained from the linear mode calculation. Note that in each 
panel of Fig. [2] we have normalized the different sets of data 
so that they have the same scale. The normalization factor 
has also been included in figure labels for easy recovery of 
the absolute value in the simuation. We can see that wave- 
functions obtained from two studies agree quite well. The 
only obvious differences are near the outer disc boundary, 
which can be attributed to the fact that outer boundary 
conditions employed in the numerical simulation and lin- 
ear mode calculation are not exactly the same. Neverthe- 



less, this agreement and the agreement in the mode growth 
rate in Fig. [1] confirm that our non- linear simulations indeed 
capture the same unstable disc p-modes as those revealed in 
linear perturbation analysis. 

To explore the time variability of the flow, we carry out 
Fourier transform of the radial velocity Ur{r,(l),t) at fixed 
r and during different evolutionary stages. In Fig. [3] we 
show some examples of the resulting power density spectra 
(normalized to the maximum value of unity). Different rows 
correspond to different total sampling times and different 
columns represent runs with different initial surface density 
perturbation. In the left columns, the disk has an initial 
perturbation with azimuthal mode number m = 3. A mix- 
ture of various modes/oscillations are excited in the flow. 
After about 10 orbits of evolution, one of them (the fastest 
growing mode) has its oscillation amplitude grown by a large 
amount such that it dominates over other modes. This corre- 
sponds to the primary spike in the top-left panel. The other 
spikes in the same panel are harmonics of this primary spike 
(see the labels of the dashed vertical lines). Table [1] shows 
that the frequenc}|f| of this fastest growing mode (cjri) dif- 
fers from the frequency obtained in linear mode calculation 
by only 0.3%, which again demonstrates the consistency of 
these two studies. After the perturbation saturates (bottom- 
left panel) , we see that the basic structure of the power den- 
sity spectrum does not change much except that these spikes 
are not as "clean" as in the linear regime; this is probably 
due to the interaction of different modes. Compared with 
the upper panel, the location of the primary spike {ujr2 in 
Table [1]) is increased by 2.0%. In Table [1] we also include 
the results from both linear mode calculation and numerical 
simulation for modes with other mode number m. The com- 
parison illustrates two main points: First, the frequencies of 
the fastest growing modes during the exponentially growth 
stage of numerical simulations are exceptionally close to the 
linear calculation results (differ by less than 1%); second, 
the frequencies of the fastest growing modes during the sat- 
uration stage are only slightly higher (except for the m = 8 
mode which shows lower frequency) than ones during the 
exponential growing stage. These indicate that the mode 
frequencies obtained in linear mode calculation are fairly 
robust and can be reliably applied in the interpretation of 
HFQPOs. 

In the right columns of Fig. [S] the simulation starts 
with a random initial perturbation which excites modes with 
various m's. During the exponential growth stage, six modes 
stand out. By examining the location of the corresponding 
spikes, we know that these are the m = 3, 4, 5, 6, 7, 8 
modes. Although the m = 3 mode seems to be the most 
prominent one in this figure, we note that this is because for 
this particular run the initial random perturbation happens 
to contain more m = 3 wave components. If we were to start 
the run with a different initial random perturbation, then 
the relative strengths of those peaks would also be different 
(not necessarily with m = 3 being dominant). 

Fig. U shows the comparison of the velocity perturba- 
tions during the linear stage (T = 20 orbits), at the end 
of the linear stage (T = 30 orbits) and during the sat- 



^ All the frequencies in this paper are angular frequencies unless 
otherwise noted. 
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Table 1. Comparison of results from linear and nonlinear studies of overstable disc p-modes 













\uJr — CJrl /^r* 


UJr2^ 


\uJr2 - ^rl|/^r2^ 


2 


1.4066 


0.0632 


0.7033 


1.3998 


0.5% 


1.4296 


2.1% 


3 


2.1942 


0.0733 


0.7314 


2.1997 


0.3% 


2.2445 


2.0% 


4 


3.0051 


0.0763 


0.7512 


2.9996 


0.2% 


3.0594 


2.0% 


5 


3.8294 


0.0751 


0.7659 


3.7995 


0.8% 


3.8886 


2.3% 


6 


4.6621 


0.0714 


0.7770 


4.6494 


0.3% 


4.7749 


2.6% 


7 


5.5007 


0.0664 


0.7858 


5.4992 


0.03% 


5.6756 


3.1% 


8 


6.3436 


0.0607 


0.7930 


6.3492 


0.09% 


6.3189 


0.5% 



^ Azimuthal mode number 

^ Mode frequency from the linear calculation (in units of Keplerian orbital frequency at the 

inner disc boundary; same for uji, cjri and u:r2) 
^ Mode growth rate from the linear calculation 

^ Ratio of wave pattern speed to the Keplerian orbital frequency at the inner disc boundary 
^ Mode frequency during the exponential growth stage of nonlinear simulation (peak frequency 

of the power density spectrum between ~ 10 orbits and ~ 30 orbits) 
^ Difference between cur (linear result) and cjri (nonlinear result) 

s Mode frequency during the saturation stage of nonlinear simulation (peak frequency of the 

power density spectrum between ~ 30 orbits and ~ 100 orbits) 
^ Difference between ujri and ujr2 
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Figure 3. Power density spectra of the radial velocity perturbations near the disk inner boundary (r = 1.1). Each panel shows the 
normalized FFT magnitude as a function of frequency. The left and right columns are for runs with initial m = 3 and with random 
perturbation, respectively. In the top and bottom panels, the Fourier transforms are sampled for time periods of [10, 30] orbits and [30, 
100] orbits, respectively. 



uration stage (T = 50 orbits) for simulations with differ- 
ent initial perturbations. At T = 20 orbits, the oscillation 
mainly comes from the single fastest growing mode and has 
a smooth radial profile. At T = 30 orbits, the perturbation 
starts to saturate, and the oscillation now consists of many 
different modes, and its radial profile exhibits sharp varia- 



tions at several locations. These sharp features remain after 
the saturation (T = 50 orbits) . 

To see this evolution from a different perspective, in 
Fig. [5] we show the color contours of radial velocity for runs 
with different initial perturbations (different columns) at dif- 
ferent times (different rows). We can clearly see that spiral 
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Figure 4. Evolution of the radial profiles of velocity perturbations from the non-linear simulations. The top and bottom panels show 
the radial and azimuthal components of velocity perturbation, respectively. The left and right panels are for cases with azimuthal mode 
number m = 2 and m = 3, respectively. In each panel, the data are taken from points with fixed cf) = 0.47r. Different line types represent 
different times during the simulation. 



waves gradually develop due to the instability. As the sys- 
tem evolves, sharp features in velcoities emerge. At the end 
(the bottom row) the spiral arms becomes more irregular, 
which is related to the emergence and interaction of multiple 
modes after saturation. 

The sharp velocity variations shown in Figs. |4]|5] suggest 
shock-like features. Note that since our simulations adopt 
isothermal equation of state (with P/S ^constant), there 
is no entropy generation and shock formation in the strict 
sense. The radial velocity jump (see Fig. |4]) is comparable 
but always smaller than the sound speed. Nevertheless, these 
sharp features imply that wave steepening plays an impor- 
tant role in the mode saturation. 



4 CONCLUSIONS 

We have carried out high-resolution, two-dimensional hydro- 
dynamical simulations of overtstable inertial-acoustic modes 
(p-modes) in BH accretion discs with various initial condi- 
tions. The evolution of disc p-modes exhibits two stages. 
During the first (linear) stage, the oscillation amplitude 
grows exponentially. In the cases with a specific azimuthal 
mode number (m), the mode frequency, growth rate and 
wavefunctions agree well with those obtained in our previous 
linear mode analysis. In the cases with random initial per- 
turbation, the disc power-density spectrum exhibits several 
prominent frequencies, consistent with those of the fastest 



growing linear modes (with various m's). These comparisons 
with the linear theory confirm the physics of corotational in- 
stability that drives disc p-modes presented in our previous 
studies (Lai & Tsang 2009; Tsang & Lai 2009c; Fu & Lai 
2011, 2012). In the second stage, the mode growth satu- 
rates and the disc oscillation amplitude remains at roughly 
a constant level. In general, we find that the primary disc 
oscillation frequency (in the cases with specific initial m) is 
larger than the linear mode frequency by less than 4%, in- 
dicating the robustness of disc oscillation frequency in the 
non-linear regime. Based on the sharp, shock-like features 
of fluid velocity profiles, we suggest that the nonlinear satu- 
ration of disc oscillations is caused by wave steepening and 
mode-mode interactions. 

As noted in Section 1, our 2D hydrodynamical simula- 
tions presented in this paper do not capture various com- 
plexities (e.g., magnetic field, turbulence, radiation) associ- 
ated with real BH accretion discs. Nevertheless, they demon- 
strate that under appropriate conditions, disc p-modes can 
grow to nonlinear amplitudes with well-defined frequencies 
that are similar to the linear mode frequencies. A number 
of issues must be addressed before we can apply our the- 
ory to the interpretation of HFQPOs. First, magnetic fields 
may play an important role in the disc oscillations. Indeed, 
we have shown in previous linear calculations (Fu & Lai 
2011,2012) that the growth of p-modes can be significantly 
affected by disc toroidal magnetic fields. A strong, large- 
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radial velocity perturbation at 20 turns 




radial velocity perturbation at 20 turns 




radial velocity perturbation at 20 turns 




Figure 5. Evolution of the radial velocity for runs with initial m = 2 (left), m = 3 (middle) and random (right) perturbations, 
respectively. From top to bottom, the times are T = 20, 30 and 50 orbits, respectively. Note that the color scale varies from panel to 
panel. 



scale poloidal field can also change the linear mode frequency 
(Yu & Lai 2013). Whether or not these remain true in the 
non-linear regime is currently unclear. Second, understand- 
ing the nature of the inner disc boundary is crucial. Our 
calculations rely on the assumption that the inner disc edge 
is reflective to incoming spiral waves. In the standard disc 
model with zero-torque inner boundary condition, the radial 
inflow velocity is not negligible near the ISCO, and the flow 
goes through a transonic point. While the steep density and 
velocity gradients at the ISCO give rise to partial wave re- 
flection (Lai & Tsang 2009), such radial inflow can lead to 
significant mode damping such that the net growth rates of 
p-modes become negative. This may explain the absence of 
HFQPOs in the thermal state of BH x-ray binaries. However, 
it is possible that the inner-most region of BH accretion discs 
accumulates significant magnetic flux and forms a magneto- 
sphere. The disc-the magnetosphere boundary will be highly 
reflective, leading to the growth of disc oscillations (Fu & 
Lai 2012; see also Tsang & Lai 2009b). Finally, the effect of 
MRI-driven disc turbulence on the p-modes requires further 
understanding. In particular, turbulent viscosity may lead 



to mode growth or damping, depending on the magnitude 
and the density-dependence of the viscosity (R. Miranda & 
D. Lai 2013, in prep). 
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